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1.  Introduction 


Recently,  Syn  and  Chen  (i)  utilized  a  small  fracture  specimen  to  study  the  fracture  behavior  of 
an  aluminum/epoxy  interface  under  high  rate  loading.  The  specimen,  due  to  its  novel  shape,  is 
frequently  referred  to  as  a  “butterfly”  specimen,  as  it  mimics  a  compact  four-point  bend 
specimen,  with  suitable  dimensions  for  high  strain-rate  testing  in  a  Kolsky  bar.  Whittie  et  al.  (2) 
utilized  a  butterfly  specimen  to  examine  the  fracture  toughness  of  cross-linked  epoxies  as  a 
function  of  loading  rate;  the  small  specimen  allowed  both  servo-hydraulic  and  Kolsky  bar 
testing.  The  specimen  geometry  was  further  developed  by  Weerasooriya  et  al.  (3)  to  study 
fracture  behavior  of  adhesive  bonds  at  low  and  high  rates  of  loading.  The  combined  research 
(1-3)  allowed  determination  of  fracture  energies  by  examining  the  load-displacement  curve  from 
either  the  servo-hydraulic  tester  or  through  traditional  split  Hopkinson  pressure  bar  (SHPB)  bar 
analysis.  However,  measurement  of  fracture  toughness  through  crack-tip  opening  displacements 
or  local  strain  fields  was  well  outside  of  the  digital  image  correlation  (DIC)  system’s  resolution 
(2,  3),  as  the  epoxies  and  adhesives  considered  were  quite  brittle  and  displacements  near  the 
crack-tip  are  on  the  order  of  10  pm.  Stress  intensity  factor  solutions  for  butterfly  specimens  are 
not  available  in  the  literature,  preventing  comparison  of  the  measured  fracture  energies  with  the 
expected  stress  intensity  factors  at  the  crack-tip. 

In  this  work,  we  sought  to  derive  a  stress  intensity  factor  solution  for  a  generalized  version  of  the 
U.S.  Army  Research  Laboratory  (ARL)  butterfly  specimen  within  the  useful  range  of  geometries 
for  Kolsky  bar  use.  The  linear  elastic  stress  intensity  factor  solution  that  has  been  investigated  in 
this  study  is  the  starting  point  to  allow  more  advanced  analysis  of  viscoelastic,  strain-rate,  and 
temperature  effects  that  are  critically  important  for  many  polymeric  materials. 


2.  Methodology 


The  finite  element  method  is  widely  used  to  calculate  stress  intensity  factors  for  various 
geometries  that  cannot  be  computed  analytically.  We  utilized  Sandia  National  Laboratories’ 
Sierra  suite  of  finite  element  codes  (4)  to  perform  quasi-static,  implicit  solid  dynamics,  and 
explicit  solid  dynamics  analyses.  The  use  of  Sierra  and  Department  of  Defense  high- 
performance  computing  assets  allowed  the  completion  of  a  multitude  of  fully  three-dimensional 
(3-D)  simulations  to  explore  the  effects  of  specimen  geometry,  loading  magnitudes  and  rates, 
frictional  effects,  and  numerical  resolution  on  the  stress  intensity  factor. 
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Computation  of  the  stress  intensity  factor  from  finite  element  results  was  performed  using  the 
J-integral  (5).  While  the  J-integral  has  documented  path-dependence  for  elastic-plastic  materials 
(particularly  in  the  plastic  zone)  (<5),  the  integral  is  path-independent  for  the  linear  and  nonlinear 
elastic  materials  considered  here.  This  path  independence  allows  the  J-integral  to  be  a  robust 
calculation  of  the  stress  intensity  factor,  as  it  does  not  rely  on  the  numerical  solution  accurately 
capturing  the  singularity  at  the  crack-tip.  The  J-integral  is  defined  in  two  dimensions  (2-D)  as 

J  =  §  (Whj  -  ajin  ju[  t  )dT  ,  ( 1 ) 

r 


where  T  represents  a  closed  path  around  the  crack-tip,  n  is  the  outward  normal  to  T,  W  is  the 
strain  energy  density,  cris  the  Cauchy  stress,  and  u  is  the  displacement.  Note  that  in  equation  1, 
the  1 -direction  is  taken  to  be  along  the  path  of  the  crack.  To  evaluate  J,  the  domain  integral 
method  (5)  is  utilized: 

JQ  =  -  j  (Wq  i  -  vfiq,  jUitl  )dA.  (2) 

A 


In  equation  2,  A  is  a  domain  that  contains  both  crack  faces  as  part  of  the  boundary,  and  the 
function  q  is  defined  such  that  q  =  1  on  the  inner  boundary  of  A  and  q  =  0  on  the  outer  boundary. 
Li  et  al.  (5)  demonstrated  that  J  =  Jq  for  nonlinear  elastic  materials;  for  inelastic  materials,  Jq  is 
an  average  value  of  J  over  paths  determined  by  curves  of  constant  q  in  the  domain  (7).  This 
domain  integral  can  be  evaluated  using  standard  finite  element  techniques.  In  this  work,  the 
function  q  was  taken  from  Carka  and  Landis  (7): 


q  = 


R„ 


R  -R 


(3) 


where  r  is  a  radial  coordinate  centered  on  the  crack-tip  and  R,  and  R„  represent  the  inner  and 
outer  radii  of  the  annular  domain  A. 


To  test  the  efficacy  of  the  numerical  simulation,  an  edge-cracked  plate  loaded  in  simple  tension 
was  simulated  in  Sierra  using  a  plane  strain  assumption.  The  plate,  with  in-plane  dimensions  of 
H  =  20  mm  and  W=  10  mm,  was  simulated  with  edge  cracks  of  lengths  a  =  1,2,  and  3  mm. 
Meshes  of  uniform  quadrilateral  elements  of  side  lengths  a/ 5,  r//15,  and  a/ 45  were  utilized.  The 
J-integral  was  evaluated,  as  described  previously,  using  the  domain  integral  for  a  variety  of  radii 
ranging  from  a/ 10  to  2a,  and  all  deviated  less  than  0.5%,  thus  verifying  path  independence. 
Recalling  that,  for  a  linear  elastic  material  under  plane  strain  deformations,  the  stress  intensity 
factor  is  related  to  the  J-integral  by 

(4) 
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one  can  compare  the  J-integral  results  to  published  values  for  the  stress  intensity  factor.  Tada  et 
al.  (5)  give  an  empirical  form  for  the  stress  intensity  factor  of  an  edge-cracked  plate  under  simple 
tension  (to  a  claimed  0.5%  accuracy)  of 


Kj  -  ct 4m (1.122-0. 231 


M 

+  10.556 

2 

-21.710 

M 

3 

+  30.382 

UJ 

UJ 

UJ 

UJ 

(5) 


with  a  representing  the  crack  length  and  b  representing  the  plate  width.  For  the  scenario  with 
a  =  2  mm  and  b  =  10  mm,  the  finite -element  results  for  the  three  mesh  resolutions  are  given  in 
table  1.  Note  the  excellent  agreement  with  the  published  results  for  the  al  15  and  a/45  cases. 
Further  meshing  trials  indicated  that  this  mesh  resolution  needs  only  to  be  maintained  in  a  small 
area  around  the  crack-tip,  dramatically  reducing  the  number  of  elements  needed  for  simulation. 


Table  1.  Stress  intensity  factors  for  various  mesh  resolutions. 


Mesh  Resolution 

Calculated  K,/(cff(-m)  °'5) 

Error 

(%) 

a/5 

1.445 

5.29 

a/ 15 

1.367 

-0.44 

a! 45 

1.360 

-0.93 

In  addition  to  the  quasi-static  plane  strain  case,  fully  3-D  dynamic  simulations  were  conducted  to 
determine  the  efficacy  of  using  the  explicit  solver  to  find  stress  intensity  factors.  The  dynamic 
effects  were  analyzed  by  evaluating  the  dynamic  J-integral  (adding  the  kinetic  energy  to  the 
strain  energy  term  in  equation  1)  as  a  function  of  time.  Monitoring  the  individual  terms  allows 
evaluation  of  the  inertial  contributions  and  time  fluctuations.  For  ramp  loads  or  pulse-to-plateau 
loadings  of  sufficiently  slow  rates,  the  contribution  of  the  dynamic  term  is  negligible  and  the  low 
strain-rate  simulation  results  are  within  numerical  error  of  the  quasi-static  analytical  results.  The 
effects  of  a  3-D  specimen  are  subtle  for  the  problem  studied  here;  the  thickness  of  the  specimen 
will  determine  if  the  results  are  closer  to  plane  stress  or  plane  strain.  Giner  et  al.  (9)  express  a 
path-area  independent  integral  Jxi  as 


J xl  J P  d  A 


^  duj^ r  d 

Wn,-<7 .. — -ft,.  dT  -  [  — 
1  "  '  J  Pbr 


dx. 


J 


A(D 


r  dut  V 

%3+-  dA- 
dxj 


(6) 


Here,  Jp  is  the  standard  line  integral  used  in  the  2-D  J-integral  with  path  T,  and  Ja  is  an  integral 
of  out-of-plane  stresses  over  the  area  A(  V)  contained  by  T.  Note  that  Jxl  is  evaluated  for  a 
position  along  the  crack-tip.  Giner  et  al.  (9)  note  that  JP  and  JA  are  a  function  of  the  position 
relative  to  the  boundaries  and  that  JA  vanishes  for  infinitesimal  T.  Jp  can  be  evaluated  using  the 
domain  integral  method  (equation  2).  JA  is  evaluated  as  an  area  integral  inside  the  surface 
bounded  by  q  =  0,  with  the  integration  weighted  by  q  (thus  giving  an  average  for  R,  <r<  Ra,  for 
the  q  in  equation  3).  For  the  edge-cracked  plate  with  H  =W=  10  mm  and  a  =  2  mm,  JxI  was 
calculated  along  the  thickness  direction  for  specimen  thicknesses  of  1,  5,  10,  20,  and  50  mm; 
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results  are  plotted  normalized  by  the  thickness  D  in  figure  1.  The  centerline  value  of  7v/ 
transitions  from  values  close  to  the  plane  stress  value  of  /  at  D  =  1  and  5  mm  to  values  close  to 
the  plane  strain  value  at  D  =  20  mm.  Also  note  that  there  are  distinct  differences  between  the 
edges  of  specimens  and  the  centerline,  which  should  be  considered  when  conducting  analyses  on 
surface-based  measurements.  That  said,  utilizing  the  explicit  solver  on  a  3-D  problem,  with 
sufficiently  slow  loading  and  long  analysis  times,  reproduces  theoretical  results  for  two 
dimensions.  The  ability  to  utilize  the  explicit  solver  allows  a  common  numerical  framework  to 
be  used  for  the  quasi-static  testing  in  this  report  to  the  planned  high  strain-rate  simulations, 
similar  to  prior  experiments  (1-3). 


x/D 


Figure  1.  Jxl  as  a  function  of  axial-coordinate  for  edge-cracked  specimens  of  thicknesses  1,  5,  10,  20, 
and  50  mm. 
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3.  ARL  Butterfly  Specimen 


While  the  experimental  measurement  of  displacement  and  force  gives  overall  fracture  energy,  it 
does  so  in  a  specific  configuration  that  is  not  readily  generalized.  Prior  research  into  similar 
geometries  (1)  utilizes  the  stress  intensity  factor  for  a  four-point  bending  specimen  without 
testing  its  efficacy  for  this  particular  nonstandard  geometry.  We  sought  to  find  a  more  accurate 
representation  of  the  stress  intensity  factor,  parameterized  for  a  generalized  version  of  this 
geometry,  utilizing  the  finite  element  method. 

Shown  in  figure  2,  the  generalized  version  of  this  geometry  has  the  following  parameters:  an 
overall  height  H  and  thickness  D.  normal  distance  from  the  centerline  to  the  upper  contact  point 
Li,  normal  distance  from  the  centerline  to  the  upper  contact  point  L2,  crack  root  distance  ao,  and 
overall  crack  length  a.  Note  that  the  radii  of  the  upper  and  lower  contact  surfaces  are  equal  to  Lj 
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3.1  Numerical  Setup  for  Butterfly  Specimen 

The  butterfly  specimen  was  developed  as  a  simple  way  to  produce  a  mode-I  crack  utilizing 
displacement-controlled  compression  tests,  particularly  for  high  strain-rate  loading.  Here,  the 
effort  was  to  characterize  the  fracture  toughness  of  a  brittle  epoxy  with  a  stiffness  of  2.2  GPa  and 
Poisson’s  ratio  of  0.3.  Numerically,  we  examined  the  compression  of  a  20-mm-deep  ( D ) 
specimen  with  metallic  platens  25  mm  in  diameter  and  6.35  mm  thick.  Recognizable  factors  that 
are  not  based  on  the  specimen  geometry  and  displacement  amplitude  arise  in  the  numerical 
simulations:  mesh  resolution,  loading  history,  platen  material,  and  frictional  effects.  To 
ascertain  the  effects  of  these  factors,  a  standard  geometry  was  selected  (L/  =  1 .7  mm; 

L2  =  4.15  mm;  H  =  13.49  mm;  ao  =  2.97  mm;  a  =  4.1  mm;  crack  is  0.3  mm  wide  with  a  rounded 
nose)  and  a  set  of  parametric  studies  was  conducted. 

Initial  simulations  were  conducted  with  steel  (E  =  200  GPa;  v  =  0.29)  platens  with  a  Coulomb 
friction  coefficient  of  0.001  on  the  standard  geometry.  It  was  found  that  a  small  friction 
coefficient  was  necessary  to  prevent  small  asymmetries  from  causing  significant  slip.  The 
specimen  was  meshed  with  a  near-uniform  mesh  resolution,  then  refined  near  the  crack-tip  and 
the  contact  areas  on  both  the  specimen  and  platen  (in  the  form  of  a  3:1  reduction  in  element  side 
length),  as  shown  in  figure  3.  It  was  found  that  the  contact  patches  needed  a  certain  level  of 
refinement  to  prevent  slip-stick  issues  for  larger  values  of  the  friction  coefficients.  To  ascertain 
the  effect  of  loading  amplitude  on  the  result,  the  standard  specimen  meshed  with  a  resolution  of 
0.6  mm  was  subjected  to  a  displacement-controlled  ramp  loading  from  0  to  1.0  mm  of 
displacement  in  10  ms.  Computation  of  Jxj  on  the  center-plane  of  the  specimen  was  conducted  at 
various  times  in  the  simulation,  and  the  load  applied  to  the  platens  was  computed  by  integrating 
the  normal  stresses  over  the  platen  surfaces.  It  was  found  that  the  applied  load  was  proportional 
to  the  square  root  of  Jxi ;  or,  in  other  words,  the  load  is  proportional  to  the  stress  intensity  factor, 
as  shown  in  figure  4. 
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Figure  3.  Sample  mesh  for  mesh  resolution  studies. 
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_ platen  load  (N) _ 

Figure  4.  Square  root  of  JxI  vs.  applied  load  for  the  standard  butterfly  specimen. 

A  mesh  resolution  study  was  conducted  to  insure  convergence  in  Jx/  before  pursuing  further 
parametric  studies.  The  lower  platen  was  held  fixed  in  the  loading  direction,  with  the  upper 
platen  given  a  constant  velocity  of  0.1  m/s.  Thus,  the  standard  mesh  shown  in  figure  3  was 
refined  by  splitting  each  node-to-node  segment  in  half  (thus  producing  double  the  mesh  density). 
Results  are  shown  in  figure  5.  There  is  less  than  a  2%  difference  between  the  Jxi  profiles  at  a 
given  platen  displacement,  echoing  the  results  from  the  edge-cracked  tension  specimen  in  section 
3.1  (and  an  attribute  of  the  nonlocal  nature  of  the  J-integral).  This  differential  is  even  smaller 
when  the  results  are  normalized  for  a  common  load  or  tip  displacement.  Further  analyses  in  this 
report  are  conducted  with  a  mesh  resolution  of  0.6  mm  (overall)/0.2  mm  (near  crack-tip). 
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x/D 


Figure  5.  Jx//2  vs.  axial  coordinate  for  a  platen  displacement  of  0.3  mm  for  standard  butterfly 
specimens  with  mesh  resolutions  of  0.6  and  0.3  mm. 

Experimentally,  aluminum  platens  are  likely  to  be  utilized  in  a  SHPB  setup  due  to  their  lower 
impedance,  which  is  important  to  obtain  a  high  signal  strength  in  the  dynamic  experiments. 
Simulations  substituting  aluminum  (E  =  68.9  GPa;  v  =  0.33)  platens  gave  no  appreciable 
differences  from  those  of  steel  platens. 

3.2  Friction  and  Specimen  Thickness  Effects 

The  effects  of  friction  were  analyzed  by  loading  the  standard  specimen  with  steel  platens;  the 
lower  platen  was  held  fixed  while  the  upper  platen  was  given  a  fixed  velocity  of  1  m/s.  Friction 
was  modeled  with  a  fixed  Coulomb  friction  coefficient  between  the  epoxy  and  steel  ranging 
between  0.01  and  1.  Initial  investigation  of  the  standard  specimen  revealed  a  nonlinear 
dependence  between  Jxl  at  a  given  platen  displacement  and  the  friction  coefficient.  Further 
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investigation  demonstrated  that  the  specimen  thickness  D  and  the  friction  coefficient  combined 
to  give  a  complicated  effect  on  Jxi.  Figure  6  shows  JxlI/2  on  the  center-plane  of  the  specimen  vs. 
the  friction  coefficient  for  specimen  thicknesses  of  5,  10,  20,  25,  and  40  mm,  at  a  platen 
displacement  of  0.3  mm.  Even  when  renormalizing  the  specimens  to  a  common 


x/D 


Figure  6.  Square  root  of  JxI  at  the  center-plane  vs.  friction  coefficient  for  several  specimen 
thicknesses. 

load-per-unit-length,  the  results  were  difficult  to  interpret.  Careful  examination  of  the  simulation 
data  revealed  that  the  displacement  difference  between  the  lower  and  upper  tips  of  the  specimen 
-  Si)  was  a  more  direct  indication  of  the  actual  load  state  of  the  crack-tip.  This  observation  is 
not  entirely  surprising;  the  displacement  difference  is  directly  related  to  the  bending  moment 
applied  to  the  specimen  and,  thus,  the  load  at  the  tip.  Figure  7  shows  Jxl1/2  at  the  specimen 
center-plane  as  a  function  of  the  tip  displacement  difference  for  various  platen  displacements 
(0.1  to  0.6  mm),  specimen  thicknesses  (5  to  50  mm),  and  friction  coefficients  (0.01  to  1.0).  Note 
the  excellent  linear  correlation  (R“  >  0.99)  between  the  tip  displacement  difference  and  Jxi  . 


10 


tip  displacement  difference  (mm) 


Figure  7.  Square  root  of  Jxl  at  the  center -plane  for  various  specimen  thicknesses,  friction 
coefficients,  and  platen  displacements. 

Specimen  thickness  is  an  important  consideration  when  analyzing  the  data  from  the  simulations. 
Figure  8  shows  Jxi  as  a  function  of  axial  position  in  the  specimen,  for  several  different  specimen 
thicknesses,  for  a  friction  coefficient  of  0.7.  Note  that  there  is  a  marked  difference  between  the 
thin  (e.g.,  5  mm)  and  the  thicker  (e.g.,  40  mm)  specimens;  this  is  similar  to  the  plane  stress  to 
plane  strain  transition  seen  in  the  edge-cracked  specimens  in  section  3.1.  The  overshoot  near  the 
edges  of  the  thicker  specimens  appears  to  be  a  frictional  effect — their  magnitude  is  reduced 
greatly  for  low  values  of  the  friction  coefficient.  Looking  at  the  center-plane  of  the  specimens, 
there  is  a  clear  transition  from  plane-stress  for  thinner  specimens  to  plane  strain  for  the  thicker 
ones.  Partitioning  the  data  in  figure  7  into  thin  (<20  mm),  transitional  (20-25  mm),  and  thick 
(>25  mm)  specimens  gives  figure  9,  where  one  can  observe  two  distinct  linear  correlations  with 
R  >  0.999.  For  a  given  tip  displacement  difference,  the  thick  correlation  ranges  from  4.45%  to 
5.1%  less  than  the  thin  correlation,  comparing  well  with  the  4.61%  predicted  for  the  plane  strain 
to  plane  stress  ratio. 
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Figure  8.  Jxj  vs.  axial  position  for  butterfly  specimens  of  various  thicknesses. 
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tip  displacement  difference  (mm) 


Figure  9.  Square  root  of  JxI  at  the  center-plane  for  various  specimen  thicknesses,  friction  coefficients,  and 
platen  displacements;  data  separated  into  thin  (<  20  mm),  thick  (>  25  mm),  and  transitional 
(20-25  mm)  sets. 

3.3  Butterfly  Specimen  Geometry 

In  order  to  generate  a  stress  intensity  factor  function  for  a  generalized  butterfly  specimen,  a  large 
number  of  analyses  were  conducted  utilizing  a  plane  strain  formulation  and  rigid  platens  with  a 
constant  friction  coefficient  (p  =  0.4).  These  assumptions  were  made  for  convenience;  the 
discovery  of  the  tip  displacements  as  an  effective  load  (section  3.2)  reduced  the  impact  of  such 
assumptions,  and  the  significantly  smaller  computational  load  of  a  plane  strain  problem  allowed 
a  large  number  of  analyses  to  be  conducted  to  compute  and  verify  the  geometric  factors. 
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The  stress  intensity  factor  was  assumed  to  have  the  general  form  of 


6  M-im  /  \ 

Ki  =7, — 

(1  -  a )  // 


where  M  represents  the  moment  (or  here  tip  displacement  differential)  and /(a)  and  g(. . .) 
represent  functions  of  nondimensional  parameters  that  characterize  the  specimen  geometry.  This 
functional  form  is  taken  from  commonly  published  forms  for  the  stress  intensity  factor  for  four- 
point  bending  specimens  (5). 


Four  nondimensional  parameters  were  selected  to  parameterize  the  geometry:  a  =  a/H 
represents  the  relative  crack  length,  L?/L  /  and  H/Lj  represent  geometric  factors  for  the  applied 
moment,  and  cio/H  represents  the  depth  of  the  wide  V  at  the  base  of  the  specimen.  There  are 
some  obvious  geometric  restrictions  (a  >  ao/H,  L2/L /  >  1,  and  H/L2  >  2),  but  practical 
considerations  and  initial  numerical  trials  revealed  realistic  ranges  for  the  nondimensional 
parameters.  The  crack  length  a  must  be  less  than  H  —  Lj  (geometrically),  but  in  practice  the 
parameter  a  should  be  in  the  range  0.2  <  a  <  0.5;  values  outside  of  this  range  see  crack-tip 
stress-fields  affected  by  the  contact  stress-fields.  The  spread  ratio  L2/L/  and  height  ratio  H/L2 
have  an  interdependent  relationship  for  practical  bounds.  Small  spread  ratios  and  large  height 
ratios  result  in  a  near-uniaxial  compression  test  rather  than  the  desired  four-point-bending-like 
test;  practically,  L2/L2  must  be  greater  than  1.75  and  H/L2  must  be  less  than  12.  The  ratio  (L2- 
L2)/H  is  a  better  indicator  of  this  effect.  Bounding  this  ratio  by  0.1  <  (L2-Lj)/H  <  0.5  prevents 
both  compression  effects  and  the  specimen  lateral  surfaces  from  affecting  the  crack-tip  stress 
field.  The  base  V  ratio  ao/H  was  found  to  have  little  effect  on  results  as  long  as  it  was 
appreciably  less  than  a  ( cio/H  <  0.7a  gave  good  results);  values  of  ao/H  approaching  a 
unsurprisingly  gave  results  resembling  a  V-notch  rather  than  a  sharp  crack. 


Initial  simulations  were  conducted  with  a  base  geometry  of  Lj  =  2  mm,  L2/L /  =  3,  H/L/  =  8,  and 
ao/H  =  2/3*a,  with  variable  a  to  ascertain  the  effects  of  a  on  the  stress  intensity  factor.  Figure 
10  shows  the  simulation  results  of  a  sweep  of  0.2  <  a  <  0.4,  normalized  by  the  tip  displacement 
difference  and  the  factor  in  equation  7.  The  fit  shown  in  figure  10  is  of  the  form 

f(cc)  =  CeA(1-a\  (8) 


where,  for  this  case,  C  =  0.005912  and  A  =  7.3045.  The  accuracy  of  the  fit  is  excellent,  with 
R2  >  0.999. 
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Figure  10.  Normalized  stress  intensity  factor  vs.  a  for  L2/Lt  =  3  and  H/Li  =  8. 


To  explore  the  other  geometric  effects,  a  was  held  constant  at  0.3  and  L2/Li  and  H/Lj  were 
varied  over  the  ranges  2  <  L2/Lj  <  3.5  and  6  <  H/Lj  <  12.  Unlike  previous  studies,  where  tip 
displacement  difference  was  the  best  indicator  of  loading,  it  was  found  that  the  fit  for  the  stress 
intensity  factor  was  more  accurate  with  a  moment  description  of 


M  oc 


S2L2  -8xL 

H 


1 


(9) 


where  S2  and  Si  are  the  displacements  of  the  lower  and  upper  tip,  respectively.  Note  that  for  a 
constant  L2/Li,  this  simply  reduces  to  the  displacement  difference  and  all  previous  results  hold. 
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Further  experimentation  with  parametric  functions  found  that  a  single  function  of  (L2-Lj)/H 
yielded  acceptable  accuracy,  as  shown  in  figure  11.  The  functional  fit  chosen  is  of  the  form 
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'l.  -lA 

C(l-B  ), 
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(10) 


where  C  =  1.3173  and  B  =  1.1093.  The  fit  prediction  was  +5%  for  all  simulated  points.  Much  of 
this  error  is  likely  attributable  to  contact  patch  discretization  differences,  as  the  scatter  of  a  set  of 
analyses  with  identical  geometry,  but  a  slightly  different  local  discretization,  is  on  this  order  of 
magnitude. 


Figure  11.  Normalized  stress  intensity  factor  vs.  (L2-L/)/H. 
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Putting  equations  7,  8,  and  10  into  a  single  function  gives 


Kj=C 


„A(l-a) 


(1  -afnH 


1  -B 


hzh 

H 


W 


J) 


(11) 


where  the  moment  is  defined  as 

M  =  E'D 2  5lLl  - S'L[  .  (12) 

H 

The  E'  in  equation  12  is  E  (Young’s  modulus)  for  plane  stress  or  E/(l-V )  for  plane  strain.  The 
specimen  thickness  D  is  taken  to  be  1  for  2-D  analyses.  The  constants  in  the  equation  are 
A  =  7.3045,  B  =  1.1093,  and  C  =  2.871  x  10"8.  The  overall  error  in  this  equation  is  bounded  by 
±5%,  much  of  which  is  numerical  setup  dependent  (when  comparing  with  computed  results). 

3.4  Crack  Width  and  Shape  Effects 

The  crack  width  and  crack  tip  shape  (from  a  to  ao)  were  held  constant  with  a  semicircular  crack- 
tip  and  a  crack  width  of  wc  =  0.3  mm  for  all  of  the  analyses  carried  out  previously.  The  a  =  0.3, 
L2/L1  =  3,  H/Li  =  8  geometry  was  simulated  with  several  crack  configurations:  (1)  crack  width 
of  0.3  mm  with  a  round  crack-tip,  (2)  crack  width  of  0.15  mm  with  a  round  crack-tip,  and  (3)  a 
perfect  “crack”  with  zero  width  and  a  sharp  crack-tip.  The  results  in  these  three  cases  varied  by 
a  maximum  of  2.1%  (the  comparison  of  the  sharp  crack  with  the  widest  notch),  indicating  that 
this  effect  is  small  for  reasonable  values  of  crack  width.  Similarly,  for  a  constant  width  notch 
with  V,  semicircular,  and  square  notch  tips,  the  difference  was  1.3%,  indicating  that  measuring 
the  exact  depth  of  the  notch  is  more  critical  than  the  exact  notch  shape. 

3.5  Validation 

To  validate  the  model,  a  random  number  generator  was  utilized  to  generate  four  random  numbers 
to  describe  the  geometry  within  the  bounds  1  mm  <  Lj  <  4  mm,  0.2  <  a  <  0.4,  2  <  L//L2  <  3.5, 
and  6  <  H/Li  <  12.  Crack  widths  were  held  to  be  a  constant  0.3  mm  with  a  semicircular  crack- 
tip.  Simulations  were  conducted  for  each  of  the  geometries  and  compared  with  the  analytical 
model.  Results  are  shown  in  table  2.  Note  that  all  errors  are  less  than  7%,  with  three  of  four 
cases  less  than  3.5%.  The  small  a  of  case  2  likely  causes  the  larger  error;  this  value  is  at  the 
extremum  of  the  calibrated  function  for  /(a). 

Table  2.  Comparison  of  analytical  and  numerical  result  for  validation  geometries. 


Case  No. 

E 

(mm) 

a 

l2/l, 

H/Li 

Jxim 

(simulation) 

((N.m)  m) 

(analytical) 
((N.m) 1/2) 

Error 

(%) 

1 

2.982 

0.2377 

2.3782 

7.1056 

21.42 

20.79 

-2.94 

2 

3.498 

0.2015 

2.6785 

9.0660 

12.57 

11.73 

-6.69 

3 

3.020 

0.2734 

2.8663 

6.1108 

27.79 

28.71 

3.33 

4 

2.215 

0.3340 

2.0311 

8.1312 

16.45 

16.26 

-1.11 
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4.  Discussion 


The  stress  intensity  factor  function  given  by  equation  1 1  gives  a  reasonably  accurate  prediction 
for  a  linear  elastic,  brittle  material  under  idealized  (2-D)  loading.  The  use  of  contact  patch 
displacements  allows  variability  in  friction  coefficients  to  be  accounted  for  in  a  straightforward 
way  with  established  DIC  techniques  that  are  already  a  part  of  the  instrumentation  suite  (2).  For 
brittle  materials,  such  as  the  epoxies  in  Whittie  et  al.  (2),  equation  11  should  provide  a 
reasonable  estimate  of  the  stress  intensity  factor  at  the  crack-tip  for  quasi-statically  loaded 
butterfly  specimens.  One  possible  use  of  this  work  is  to  design  specimen  geometries  tailored  to 
the  properties  of  new  epoxies  (stiffness,  strength,  fracture  toughness,  etc.). 

The  analysis  in  this  work  has  some  limitations.  The  linear  elastic  constitutive  relation  is,  of 
course,  a  simplification.  Obvious  ramifications  are  the  lack  of  rate  effects  (strain-rate  hardening, 
viscoelasticity,  and  fracture  energy  as  a  function  of  rate),  but  it  also  eliminates  local  yielding 
from  the  analysis.  This  fact  should  be  considered  when  designing  specimens,  as  the  local 
stresses  from  the  small  contact  patches  can  be  quite  large,  leading  to  yielding  at  the  contact 
points  before  fracture  begins.  Yielding  can  also  occur  at  the  crack-tip  before  fracture  begins, 
although  typical  plastic-zone  techniques  developed  in  the  fracture  mechanics  community  can 
compensate  for  this  yielding.  Note  that  while  the  finite  element  analysis  includes  inertia,  the 
loading  pulses  are  long  and  well  conditioned.  When  conducting  dynamic  analyses,  care  must  be 
taken  to  give  the  specimen  mono  tonic  loading  (e.g.,  through  the  use  of  pulse  shaping  on  the 
SHPB),  as  vibrational  modes  of  the  specimen  can  be  excited  under  other  loading  types,  leading 
to  hard-to-interpret  results. 

Future  extensions  of  this  work  will  include  rate-dependent  effects  in  both  the  material 
constitutive  behavior  and  the  fracture  toughness  of  the  material,  and  application  of  similar 
analysis  for  the  use  of  the  ARL  butterfly  geometry  to  study  the  fracture  behavior  in  adhesive 
bonds.  Whittie  et  al.  (2)  noted  that  the  fracture  energy  nominally  doubled  when  the  loading  rate 
was  increased  from  7  to  70  kN/s.  Whether  this  change  in  fracture  energy  is  a  result  of  strain-rate 
hardening,  viscoelasticity,  or  rate  dependency  of  fracture  energy  remains  to  be  seen. 
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